The Pairing Interaction in the 2D Hubbard Model 
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A dynamic cluster quantum Monte Carlo approximation is used to study the effective pairing 
interaction of a 2D Hubbard model with a near neighbor hopping t and an on-site Coulomb interac- 
tion U . The effective pairing interaction is characterized in terms of the momentum and frequency 
dependence of the eigenfunction of the leading eigenvalue of the irreducible particle-particle vertex. 
The momentum dependence of this eigenfunction is found to vary as {coskx — cosfc^) over most 
of the Brillouin zone and its frequency dependence is determined by the exchange energy J. This 
implies that the effective pairing interaction is attractive for singlets formed between near-neighbor 
sites and retarded on a time scale set by . The strength of the pairing interaction measured by 
the size of the d-wave eigenvalue peaks for U of order the bandwidth 8t. It is found to increase as 
the system is underdoped. 



I. INTRODUCTION 

Results from numerical studies suggest that the two- 
dimensional (2D) Hubbard model exhibits the basic phe- 
nomena which are seen in the high- Tc cuprate materials. 
At half-filling, one finds a groundstate with long-range 
antiferromagnetic order^. Doped away from half-filling 
there is a pseudogap regime^"^, and at low temperature 
a striped phase^'^ as well as d-wave pairing-'^". In addi- 
tion, the various phases appear delicately balanced with 
respect to changes in parameters. All of these features 
remind one of the actual cuprate materials, so that it is 
of interest to understand the structure of the interaction 
that leads to c?a;2_j,2-wave pairing in this model. 

Here we will study this interaction for a 2D Hubbard 
model^^ on a square lattice which contains a one-electron 
near-neighbor hopping t and an onsite Coulomb interac- 
tion U. In this case the results depend upon just two 
parameters U/t and the average site occupation (n) . Pre- 
vious work""^^ has shown that the pairing interaction with 
U/t = 4 and (n) = 0.85 increases with momentum trans- 
fer, has a Matsubara frequency dependence similar to 
that of the Q = (tt, tt) spin susceptibility and is medi- 
ated by a particle-hole S = 1 exchange channel. Here we 
extend this work to explore the pairing interaction for 
larger values of U/t and various dopings. We are par- 
ticularly interested in the case in which U is of order 
of the bandwidth 8t. In this case, well developed upper 
and lower Hubbard bands are seen in the single particle 
density of states 
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Fig. 1 shows N{lu) with U/t — 8 for various fillings. 
These results were obtained from a maximum entropy 
continuation^'^ of dynamic cluster quantum Monte Carlo 



data on a 16 site cluster. At half-filling, Fig. la, iorTr^t 
one sees broad upper and lower Hubbard bands. Then, as 
the temperature scale drops below the exchange energy 
scale J ~ At^ /U , two additional structures appear above 
and below u = Q. As Preuss et. al^'^ first showed, the 
inner structures arise from the formation of two coherent 
bands, each of width ^ 2J that form as the antiferro- 
magnetic correlations develop. This characteristic Mott- 
Hubbard- antiferromagnetic four band structure^'^'^^ is 
clearly seen when U becomes of order or larger than the 
bandwidth. 

Fig. lb shows N{uj) for the doped (n) — 0.9 case. 
Here the chemical potential has moved down into the 
lower coherent band and the spectral weight in the up- 
per coherent band has essentially vanished. Nevertheless, 
the remnants of the upper and lower Hubbard bands re- 
main. Part of the motivation for this study is to examine 
the structure of the pairing interaction in this parameter 
regime. 

As in our previous work, the pairing interaction F^*^ 
will be calculated using a dynamic cluster quantum 
Monte Carlo simulation^^. As shown in Fig. 2, F^^ is 
the irreducible part of the 4-point particle-particle ver- 
tex in the zero center of mass and energy channel. Pre- 
viously, we discussed how one could extract r^'^'(fc|fc') = 
rPP{k, ~k; k' , —k') with k — (k, zwn) using a dynamic 
cluster quantum Monte Carlo simulation. 

While we studied F^^ in our earlier work, here we will 
focus on the momentum and Matsubara frequency depen- 
dence of the d^2_y2--wave eigenfunction <I>rf(k, Wn) of the 
homogeneous particle-particle Bethe-Salpeter equation 

E ^"'iMk') G^{k') Gx(-fc') $,(fc') = XMk) . 

fe' 

(2) 

At low temperatures, $d(k, cj„) has the largest eigen- 
value and this eigenvalue goes to one dX T = Tc- As 
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FIG. 1: The single particle density of states N{uj) versus ui 
for U = 8t, Nc = 15 and various values of the temperature T 
for site fillings of a) (n) = 1.0, b) (n> = 0.90. 
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FIG. 2: The Bethe-Salpeter equation for the particle-particle 
channel showing the relationship between the four-point ver- 
tex r and the particle-particle irreducible vertex F''''. The 
solid lines are dressed single particle Green's functions. 



T approaches Tc, the momentum and frequency depen- 
dence of <i>rf(k, LUn) reflect the structure of the pairing in- 
teraction at Tc, just as the superconducting gap function 
reflects the k and lu dependence of the pairing interac- 
tion in the superconducting state. Thus, while $ci(k,a;„) 
is not a quantity that is directly measurable, it has a k- 
dependence related to the momentum dependence of the 
interaction and a Matsubara frequency dependence which 
decays beyond a characteristic frequency associated with 
the dynamic character of the interaction. It also has the 



great advantage of depending upon one momentum and 
frequency variable as opposed to the multiple momentum 
and frequency variables of TPP{k\k'). 

In the following section II we review the dynamic 
cluster approximation and discuss how one calculates 
the c?a;2_y2-wave eigenvalue XdiT) and eigenfunction 
<I>d(k,w„). Then, in Sec. Ill we investigate how Xd{T) 
depends upon U and (n). Following this, we examine 
the k-dependence of $d(k, w„) and see how closely it fol- 
lows the simple cos(A:a;) — cos{ky) dependence. If it were 
of this form over the entire Brillouin zone, then it would 
imply a strictly near-neighbor pairing interaction. Then 
we turn to the a;„-dependence which reflects the dynam- 
ics of the pairing interaction and study its dependence on 
U. In Sec. IV, based upon the results for <I>rf(k, Wn), we 
construct a simple separable representation of rPP{k\k') 
and discuss the strength of the pairing interaction. Sec. V 
contains our conclusions. 



II. THE DYNAMICAL CLUSTER 
APPROXIMATION 

The Dynamical Cluster Approximation (DCA)^^ maps 
the bulk lattice to a finite size cluster embedded in a 
self-consistent bath designed to represent the remain- 
ing degrees of freedom. Short-range correlations within 
the cluster are treated explicitly, while the longer-ranged 
physics is described by a mean-field. By increasing the 
cluster size, the DCA systematically interpolates between 
the single-site dynamical mean-field result and the exact 
result, while remaining an approximation to the thermo- 
dynamic limit for finite cluster size. 

The essential assumption is that short-range quanti- 
ties, such as the single-particle self-energy E, and its 
functional derivatives, the two-particle irreducible vertex 
functions, are well represented as diagrams constructed 
from a coarse-grained propagator G. To define G, the 
Brillouin zone in two dimensions is divided into Nc = 
cells of size iTijl?. As illustrated in Fig. 3, each cell 
is represented by the cluster momentum K in its center. 
The coarse-grained Green function G(K) is then obtained 
from an average over the N /N^ wave- vectors k within the 
cell surrounding K, 



G(K,c^„) 
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Here the self-energy for the bulk lattice I](K-f k, a;„) has 
been approximated by the cluster self-energy Ec(K,ti;„). 
Consequently, the compact Feynman diagrams con- 
structed from G'(K,Ci;„) collapse onto those of an effec- 
tive cluster problem embedded in a host which accounts 
for the fluctuations arising from the hybridization be- 
tween the cluster and the rest of the system. The non- 
interacting part of the effective cluster action is then de- 
fined by the cluster-excluded inverse Green's function 



(K, Lo^) = G-i(K, c^„) + S,(K, oj„ 



(4) 
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FIG. 3: In the DCA the Brillouin zone is divided into Nc cells 
each represented by a cluster momentum K. Irreducible quan- 
tities such as the single-particle self-energy E and two-particle 
irreducible vertex F are constructed from coarse-grained prop- 
agators G'(K) that are averaged over the momenta k' within 
the cell represented by K. The cluster momenta K for a 4-site 
cluster are shown on the left and those for a 24-site cluster on 
the right. 



which accounts for the hybridization between the cluster 
and the host. Given Q^^(K.,LUn) and the interaction on 
the cluster U '^^rii^riii^ one can then set up a Hirsch- 
Fye quantum Monte Carlo algorithm^^ to calculate the 
cluster Green function and from it the cluster self-energy 
Ec(K,a'„) which is used in Eq. (3) to re-calculate the 
coarse-grained Green function (^(K, cj„)^^'^^. This pro- 
cess is then iterated to convergence. 

Since a determinantal Monte Carlo method is used, 
there is also a sign problem for the doped Hubbard model. 
However, the coupling of the cluster to the self-consistent 
host significantly reduce the sign problem so that lower 
temperatures can be reached"'^^. 

The DCA cluster one- and two-particle Green's func- 
tions that we calculate have the standard finite temper- 
ature definitions 



(5a) 



and 



Gc2cr4---cri{X4, X3; X2, Xl) — 

- (T,C,,(X4)c,3(^3)4(^2)4(^l)) . (5b) 

Here, Xf — (X£,T£), where denotes a site in the 
DCA cluster, T£ the imaginary time, Tr is the usual r- 
ordering operator, and cJ.(X2) creates a particle on the 
cluster with spin a. Fourier transforming on both the 
cluster space and imaginary time variables gives Gc{K) 
and Gc2iK4,K3;K2,Ki) with K = (K,itj„,cr). Us- 
ing Gc{K) and Gc2{K4, K3; K2, Ki), one can extract the 
cluster four-point vertex T from 

GMKi, K3; K2, Ki) = 

— Gc{Ki)Gc{K2) [5Ki,Ki&K2,K3 ^ ^Ki,K3^K2,Ki] 
T 

^Ki+Ki.Ks+KiGciKi) Gc{K:i)V {K4, i^a; i^2, Ki) 
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X GrXK2)GciKl). 



Then, using Gc and F, one can determine the ir- 
reducible particle-particle vertex F''*' from the Bethe- 
Salpeter equation shown in Fig. 2. 

Using FPP, the eigenvalues and eigenfunctions of the 
Bethe-Salpeter equation may then be calculated from 

T 



X G^ik')Gii-k')<j,^iK') ^ Kq^^iK) (7) 



Here, the sum over k' denotes a sum over both momen- 
tum k' and Matsubara w„' variables. We decompose 
k' = K' -f k'. By assumption, irreducible quantities like 
and (pa do not depend on k', allowing us to coarse- 
grain the Green function legs, yielding an equation that 
depends only on coarse-grained and cluster quantities 



|- ^ FPP (A^ -K; K', -K') xl^K') ^K') = 

KMK) (8) 



K' 



(6) 



with xl^iK') = ^Ek'GT(K' + k\tu;n')Gi{-K' - 

Here, we show a number of results for the 4-site cluster 
shown on the left in Fig. 3, which allows us to investigate 
larger values of U and lower temperatures than the 24- 
site cluster. As discussed in Ref.^*', the 4-site cluster does 
not allow for the effect of pairfield phase- fluctuations. 
Simulations on the 24-site cluster were used to determine 
the k dependence of <i>d(k, u;„). 



III. THE STRUCTURE OF THE PAIRING 
INTERACTION AS REFLECTED IN Xd{T) AND 
*d(k, 

The temperature dependence of the d-wave eigenvalues 
Xd{T) calculated using a 4-site cluster for a site filling 
(n) — 0.9 with U/t — 4, 8 and 12 are shown in Fig. 4a. 
Fig. 4b shows Ad versus U/t for T = 0.15t. Here one 
sees that it is favorable to have a Coulomb interaction 
strength U of order the bandwidth 8t. This is consistent 
with the notion that it is important to have strong short- 
range antiferromagnetic correlations. However, because 
the exchange interaction J ~ At^ /U at strong coupling, 
the short-range antiferromagnetic correlations decrease 
when U becomes large compared to the bandwidth. 

The dependence of \d{T) on the filling (n) is illustrated 
in Fig. 5 for U/t = 6. What one sees is that Xd{T) in- 
creases as the system is doped towards half-filling. How- 
ever, at half-filling the dominant eigenvalue of the 4-point 
vertex occurs in the Q ~ (tt, tt) irreducible particle-hole 
S = 1 channel as indicated by the dashed line in Fig. 5, 
and the groundstate at T = has long-range antiferro- 
magnetic order. 

For (n) = 0.9 and U/t — 8, the momentum dependence 
of the d-wave eigenvector $rf(K, w„) for the 24-site clus- 
ter at ujn = ttT is shown in Fig. 6. Here, for T/t = 0.22, 
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(a) (n) = 0.90; Nc = 4: 
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FIG. 4: (a) The d^2_y2 eigenvalue Xd{T) versus T/t for U = 
4t, 8t and 12t and (n) = 0.90. (b) The d^2_y2 eigenvalue 
Ad(T) versus U/t for T = 0.15t and (n> = 0.90. 



Xd{T) = 0.42 and the values of K lay along the dashed 
line shown in Fig. 3. One clearly sees the d-wave struc- 
ture of ^d- The dependence of <i>c;(K,7rT) for K along 
the Kx axis is shown in the inset of Fig. 6. Here, one 
sees that $(i(K, ttT) falls off as K moves away from the 
Fermi surface towards the zone center. 

We have also calculated the projection of <l>rf(K,7rT) 
on the first and second (1^2 _y2 crystal harmonics 



(9) 



(n> = 1.00 
(n) = 0.90 
(n) = 0.85 
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FIG. 5: The d^2_y2 eigenvalue XdiT) versus T/t for various 
band fillings (n) for U/t = 6. The dashed line represents the 
the leading eigenvalue Xaf in the Q = (tt, tt), S = 1 particle- 
hole channel at half-filling. 
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FIG. 6: The d^2_y2 eigenvector "l>d(K,ciJ„) at tj„ — nT, nor- 
malized to its value at K = (tt, 0), versus K for U/t = 8, band 
filling (n) = 0.9 and T/t = 0.22. In the main figure, the K 
points move along the dashed line shown in Fig. 3. The inset 
shows the behavior of $d when K varies along the axis. 



since 52 (K) = on the momenta K along the dashed 
line. 



K 



with (7i(K) = cosKx — cos Ky and 52 (K) 
cos2Ky. In table I, we list the values of d2/di versus 
[/ at a filling (n) — 0.9. Here the sum in Eq. (9) is 
over the entire Brillouin zone and the temperature was 
adjusted so that the d-wave eigenvalue for each U/t 
was the same (Ad w 0.4). If the sum over K in Eq. (9) 
is restricted to values which lay along the dashed line in 
Fig. 3, this ratio vanishes exactly in the 24-site cluster. 





u/t 


4 6 8 


cos2iir2; — 


d2/di 


0.064 0.128 0.157 



TABLE I: The ratio of the second to the first crystal d-wave 
harmonic projection of "I>d(K,7rr) for (n) = 0.9 and Xd ~ 0.4. 



The Matsubara frequency 
$rf(K,u;„)/$<j(K,7rT) with K = 
Fig. 7 for (n) = 0.9 and U/t = 4, 8 



dependence of 
(vr, 0) is shown in 
and 12. Also shown 
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in each case is the frequency dependence of x(Q,a;„) 
with Q = (tt, tt). The decrease of the characteristic 
exchange energy as U /t increases is seen in Figs. 7a-c. 
It is clear from these results that the dynamics of the 
pairing interaction T^p ^ reflected in the Un dependence 
of <i>(i(k;, cjn), is associated with the spin-fluctuation 
spectrum. 



IV. A SEPARABLE REPRESENTATION OF 
THE PAIRING INTERACTION 

With the results obtained for the d-wave eigenvector 
$(i(_ftr), one can construct a simple separable representa- 
tion of the pairing interaction T'P'p{K\K') 



(a) U = 4f; N, 



TPP{K\K') = ~Vd^d{K) ^d{K') 



(10) 



Using this separable form for Vpp{K\K'), one finds that 
Eq. (8) gives 



(11) 



K' 



Then, using the dressed DCA Monte Carlo single particle 
Green's functions and the DCA results for Ad, we can de- 
termine the strength Vd of the separable interaction from 
Eq. (11). The strength of Vd depends upon both the site 
occupation (n) and the temperature T. An alternative 
way of extracting a simple approximation for the pair- 
ing interaction is to use the d-wave projected irreducible 
vertex for ujr, = lu,-,' = ttT 



(r) 



^ ^ .g(K)rPP(K, ttTIK', ^T)g(K') (12) 



K,K' 



with (7(K) = {cosKx — cosKy)/2. 

Fig. 8 shows Vd{T) and vj^^ (T) for U = 8t and various 
values of the site occupation (n). Both approximations 
give very similar results. It is interesting to see that Vd 
becomes stronger as (n) approaches half-filling. This is 
characteristic of the Mott-Hubbard system and has been 
previously observed. A determinantal quantum Monte 
Carlo calculation of the d-wave eigenvalue on an 8 x 8 lat- 
tice at half-filling with U — 8t found that Ad approached 
1 as the temperature was lowered^^. As noted there, in 
this case the antiferromagnetic eigenvalue remained dom- 
inant and on an infinite lattice the groundstate would 
have long-range antiferromagnetic order at T = 0. A re- 
lated behavior is seen for the two-leg ladder, where the 
pair binding energy of a finite ladder is greatest for the 
first two holes that are added^°. This latent pairing ten- 
dency of the half-filled Hubbard model is also seen in 
the magnitude of the probability amplitude for adding 
two near-neighbor holes in a d-wave state reported by 
Plekhanov et al.'^^ . 

Although Vd increases as (n) goes to 1, the number of 
holes that are available for pairing is suppressed as is Tc- 



*d((7r,0),tu„) 
Xa ((7r,7r),aj™) 




n(m) / 



FIG. 7: The Matsubara frequency dependence of 
<E>d(K,t^„)/$d(K,7rT) with K = (7r,0) for (a) U/t = 4, 
(b) U/t = 8 and (c) U/t = 12 calculated for Nc^4. 
Here uj„ = (2n + 1)ttT with T/t = 0.125 and the 
band filling (n) — 0.90. Also shown is the fre- 
quency dependence of the normalized spin susceptibility 
2x(Q,a;^)/[x(Q,0) + x(Q,27rr)] for Q = (7r,7r). 
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FIG. 8: The pairing interaction Vd versus temperature for 
different dopings for U = 8t calculated for Nc = 4. 



0.030 r 



0.025 ■ 



0.020 



0.015 ■ 



0.010 ■ 



0.005 ■ 



0.000 




A measure of this is given by 



N, 



(13) 
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FIG. 9; Plot of PdoiT) versus T for [/ = 8t and (n) = 1 
showing the effect of the opening of the Mott-Hubbard gap 



which is plotted in Fig. 9 for U = 8t and (n) = 1. Here 
one clearly sees that as the temperature is lowered and 
the Mott-Hubbard gap opens, Pdo{T) is suppressed. 



V. CONCLUSION 

The cosfc^; — cosfcj, dependence of $d(k, ix>„) reflects 
a pairing interaction r''P(k|k') which increases at large 
momentum transfer k — k', implying a spatially short- 
range interaction which is repulsive for pair formation 
on the same site but attractive for singlet pair forma- 
tion between near-neighbor sites. The lu„ dependence of 
<i>rf(k, Un) tells us that the pairing interaction is retarded 
on a time scale set by J^^. The strength of the inter- 
action is largest when U is of order the bandwidth and 
increases as the system is doped towards half-filling. Of 
course, with U — 8t, as (n) goes to 1, a Mott-Hubbard 
gap opens and there are no holes to pair. 
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